from numpy import *
from numpy import tensordot as td

a1=linspace(0,10,3)
b1=linspace(3,6,3)
a2=arange(0.0,9.0,1.0)
b2=linspace(0.0,10.0,9)
c2=linspace(5.0,12.0,9)
d2=linspace(0.0,1.0,9)
e2=linspace(0.0,3.0,9)
f2=linspace(3.0,4.0,9)
a4=linspace(0,10,81)
b4=linspace(1,12,81)
c4=linspace(2,3,81)
a5=linspace(2,18,243)
a6=linspace(0,30,729)
s0=linspace(1,10,36)
m0=linspace(0,3,36)


a2.shape=(3,3)
b2.shape=(3,3)
c2.shape=(3,3)
d2.shape=(3,3)
e2.shape=(3,3)
f2.shape=(3,3)
a4.shape=(3,3,3,3)
b4.shape=(3,3,3,3)
c4.shape=(3,3,3,3)
a5.shape=(3,3,3,3,3)
a6.shape=(3,3,3,3,3,3)
s0.shape=(3,12)
m0.shape=(3,12)
kdel = identity(3)

######################################################################
##                        JUST OUTER EXAMPLES                       ##
######################################################################

###################################################
##  ans_{ijk}=a_{ij}b_{k}                        ##
###################################################
ffor = zeros((3,3,3))

for i in range(3):
  for j in range(3):
    for k in range(3):
      ffor[i,j,k] = a2[i,j]*a1[k]
      
fout = outer(a2,a1).reshape((3,3,3))

allclose(ffor, fout)

###################################################
##  ans_{ijkl}=a_{ij}b_{kl}                      ##
###################################################
gfor = zeros((3,3,3,3))

for i in range(3):
  for j in range(3):
    for k in range(3):
      for l in range(3):
        gfor[i,j,k,l] = a2[i,j]*b2[k,l]
      
gout = outer(a2,b2).reshape((3,3,3,3))

print allclose(gfor, gout)

###################################################
##  ans_{ijklmn}=a_{ij}b_{kl}c_{mn}              ##
###################################################
hfor = zeros((3,3,3,3,3,3))

for i in range(3):
  for j in range(3):
    for k in range(3):
      for l in range(3):
        for m in range(3):
          for n in range(3):
            hfor[i,j,k,l,m,n] = a2[i,j]*b2[k,l]*c2[m,n]
      
hout = outer(outer(a2,b2),c2).reshape((3,3,3,3,3,3))

print allclose(hfor, hout)

################################################################
##  ans_{ijklmnopqrstu}=a_{ij}b_{kl}c_{mn}d_{op}e_{qr}f_{st}  ##
################################################################
ifor=zeros(linspace(3,3,12))
for i in range(3):
  for j in range(3):
    for k in range(3):
      for l in range(3):
        for m in range(3):
          for n in range(3):
            for o in range(3):
              for p in range(3):
                for q in range(3):
                  for r in range(3):
                    for s in range(3):
                      for t in range(3):
                        ifor[i,j,k,l,m,n,o,p,q,r,s,t]=a2[i,j]*b2[k,l]*c2[m,n]*d2[o,p]*e2[q,r]*f2[s,t]
      
iout = outer(outer(outer(outer(outer(a2,b2),c2),d2),e2),f2).reshape(linspace(3,3,12))

print allclose(ifor, iout)

######################################################################
##                  OUTER WITH MULTIPLICATION                       ##
######################################################################


###################################################
##  ans_{ijkl}=C_{mnop}a_{im}b_{jn}c_{ko}d_{lp}  ##
###################################################
dfor=zeros((3,3,3,3))
#big tensordot example
for i in range(3):
  for j in range(3):
    for k in range(3):
      for l in range(3):
        for m in range(3):
          for n in range(3):
            for o in range(3):
              for p in range(3):
                dfor[i,j,k,l]+=a4[m,n,o,p]*a2[i,m]*b2[j,n]*c2[k,o]*d2[l,p]

dtemp1 = outer(outer(a2,b2),outer(c2,d2)).reshape(4*shape(a2))
dout = td(dtemp1, a4,((1,3,5,7),(0,1,2,3)))

print allclose(dfor, dout)

###################################################################
##  ans_{ijkl}=C_{opqr}a_{io}b_{jp}c_{kq}d_{lr}                  ##
###################################################################
jfor=zeros(linspace(3,3,4))

for i in range(3):
  for j in range(3):
    for k in range(3):
      for l in range(3):
        for o in range(3):
          for p in range(3):
            for q in range(3):
              for r in range(3):
                jfor[i,j,k,l] += a2[i,o]*b2[j,p]*c2[k,q]*d2[l,r]*a4[o,p,q,r]

jtemp1 = outer(outer(outer(a2,b2),c2),d2).reshape(linspace(3,3,8))
jout = td(jtemp1, a4,((1,3,5,7),(0,1,2,3)))

print allclose(jfor, jout)


###################################################################
##  ans_{ijklm}=C_{opqr}a_{io}b_{jp}c_{kq}d_{lr}e_{ms}           ##
###################################################################
kfor=zeros(linspace(3,3,5))

for i in range(3):
  for j in range(3):
    for k in range(3):
      for l in range(3):
        for m in range(3):
          for o in range(3):
            for p in range(3):
              for q in range(3):
                for r in range(3):
                  for s in range(3):
                    kfor[i,j,k,l,m] += a2[i,o]*b2[j,p]*c2[k,q]*d2[l,r]*e2[m,s]*a5[o,p,q,r,s]

ktemp = outer(outer(outer(outer(a2,b2),c2),d2),e2).reshape(linspace(3,3,10))
kout = td(ktemp, a5,((1,3,5,7,9),(0,1,2,3,4)))

#CANNOT USE == BECAUSE NUMERICAL ERROR IS MANIFESTED BY THIS MANY OPERATIONS, ERROR ON ORDER OF E-11
print allclose(kfor, kout)


###################################################################
##  ans_{ijklmn}=C_{opqrst}a_{io}b_{jp}c_{kq}d_{lr}e_{ms}f_{nt}  ##
###################################################################
efor=zeros(linspace(3,3,6))
#big tensordot example
for i in range(3):
  for j in range(3):
    for k in range(3):
      for l in range(3):
        for m in range(3):
          for n in range(3):
            for o in range(3):
              for p in range(3):
                for q in range(3):
                  for r in range(3):
                    for s in range(3):
                      for t in range(3):
                        efor[i,j,k,l,m,n]+=a6[o,p,q,r,s,t]*a2[i,o]*b2[j,p]*c2[k,q]*d2[l,r]*e2[m,s]*f2[n,t]

etemp1 = outer(outer(outer(outer(outer(a2,b2),c2),d2),e2),f2).reshape(linspace(3,3,12))
eout = td(etemp1, a6,((1,3,5,7,9,11),(0,1,2,3,4,5)))

#CANNOT USE == BECAUSE NUMERICAL ERROR IS MANIFESTED BY THIS MANY OPERATIONS, ERROR ON ORDER OF E-11
print allclose(efor, eout)


######################################################################
##            MORE SPECIFIC OUTER WITH MULTIPLICATION               ##
######################################################################
###################################################################
##  ans_{ij}=a_{ki}b_{kj}                                        ##
###################################################################
lfor=zeros(linspace(3,3,2))

for i in range(3):
  for j in range(3):
    for k in range(3):
      lfor[i,j]+=a2[k,i]*b2[k,j]

lout = td(a2,b2,(0,0))

print allclose(lfor, lout)

###################################################################
##  ans=a_{1i}b_{ij}c_{j1}                                       ##
###################################################################
mfor=0.0

for i in range(3):
  for j in range(3):
    mfor+=a2[0,i]*b2[i,j]*c2[j,0]

mtemp = outer(a2[0,:],c2[:,0]) #temp_{ij}=a_{1i}c_{j1} (Form same rank matrix for multiplication)
mout = td(mtemp, b2, ((0,1),(0,1)))

print allclose(mfor, mout)

###################################################################
##  ans_{ij}=a_{jl}b_{lm}del_{1i}c_{1m}                       ##
###################################################################
ofor=zeros(linspace(3,3,2))

for i in range(3):
  for j in range(3):
    for l in range(3):
      for m in range(3):
        ofor[i,j]+=a2[j,l]*b2[l,m]*kdel[0,i]*c2[0,m]

#step by step
otemp1 = outer(kdel[0,:],c2[0,:]) #temp1_{im}=del_{1i}c2_{1m}
otemp2 = td(otemp1,b2,(1,1)) #temp2_{il}=temp1_{im}b2_{lm}
oout = td(otemp2,a2,(1,1)) #out_{ij}=temp2_{il}a2_{jl}

#single expression
oshort=td(td(outer(kdel[0,:],c2[0,:]),b2,(1,1)),a2,(1,1))

print allclose(ofor, oout)
print allclose(ofor, oshort)

######################################################################
##          OUTER, MULTIPLICATION, AND WHOLE EXPRESSIONS            ##
######################################################################
###################################################################
##  ans_{ij}=(del_{ik}+del_{1i}del_{1k))a_{kj)                   ##
###################################################################
pfor=zeros(linspace(3,3,2))

for i in range(3):
  for j in range(3):
    for k in range(3):
      pfor[i,j]+=(kdel[i,k]+kdel[0,i]*kdel[0,k])*a2[k,j] 

#step by step
ptemp1 = outer(kdel[0,:],kdel[0,:]) #temp1_{ik}=del_{1i}del_{1k}
ptemp2 =  kdel + ptemp1 #temp2_{ik}=del_{ik}+temp1_{ik}
pout = td(ptemp2,a2,(1,0)) #out_{ij}=temp2_{ik}a2_{kj}

#single expression
pshort = td((kdel+outer(kdel[0,:],kdel[0,:])),a2,(1,0))

print allclose(pfor, pout)
print allclose(pfor, pshort)

###################################################################
####        Form of Equation 46 with no slip systems           ####
###################################################################
##  ans=a_{ig}b_{gk}a_{jk}a_{ib}a1_{b}b1_{d}c_{dj}
###################################################################
qfor=0.0

for i in range(3):
  for j in range(3):
    for k in range(3):
      for b in range(3):
        for d in range(3):
          for g in range(3):
            qfor+=a2[i,g]*b2[g,k]*a2[j,k]*a2[i,b]*a1[b]*b1[d]*c2[d,j]

#step by step
qtemp1 = td(a2,b2,(1,0)) #temp1_{ik}=a2_{ig}b2_{gk}
qtemp2 = td(qtemp1, a2, (1,1)) #temp2_{ij}=temp1_{ik}a2_{jk}
qtemp3 = td(a2, qtemp2, (0,0)) #temp3_{bj}=a2_{ib}temp1_{ij}
#qtemp3 = td(qtemp2, a2, (0,0)) #temp3_{bj}=temp2_{ij}a2_{ib} !NOT THIS ORDER, THIS PRODUCES temp3_{jb}
qtemp4 = td(qtemp3, c2, (1,1)) #temp4_{bd}=temp3_{bj}c2_{dj}
qtemp5 = outer(a1,b1) #temp5_{bd}=a1_{b}b1_{d}

qout = td(qtemp4, qtemp5, ((0,1),(0,1))) #out = temp4_{bd}temp5_{bd}

#single expression
qshort = td(td(td(a2,td(td(a2,b2,(1,0)),a2,(1,1)),(0,0)),c2,(1,1)),outer(a1,b1),((0,1),(0,1)))

print allclose(qfor, qout)
print allclose(qfor, qshort)

#############################################################################
####          Form of Equation 46 with slip systems                      ####
#############################################################################
##  ans^{alpha}=a_{ig}b_{gk}a_{jk}a_{ib}s0_{b}^{alpha}m0_{d}^{alpha}c_{dj} ##
#############################################################################
rfor=zeros(linspace(12,12,1))
for w in range(12):
  for i in range(3):
    for j in range(3):
      for k in range(3):
        for b in range(3):
          for d in range(3):
            for g in range(3):
              rfor[w]+=a2[i,g]*b2[g,k]*a2[j,k]*a2[i,b]*s0[b,w]*m0[d,w]*c2[d,j]

rshort=zeros(linspace(12,12,1))
for w in range(12):
  rshort[w] = td(td(td(a2,td(td(a2,b2,(1,0)),a2,(1,1)),(0,0)),c2,(1,1)),outer(s0[:,w],m0[:,w]),((0,1),(0,1)))

print allclose(rfor, rshort)






















